/*
 * TRISpH.cpp
 *
 *  Assembled on: 16 Aug 2011
 *      Author: Allan
 */

#include "TRISpH.h"

// C++ includes
#include <iomanip>
#include <iostream>
#include <memory>
#include <string>
#include <vector>
using namespace std;

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "GeoReact.h"
#include "PitzerTris.h"

const Multiphase
AssembleTRISMultiphase();

const ReactionSystem
AssembleTRISReactionSystem(const Multiphase& multiphase);

const GeoSystem
AssembleSystemTRIS();

void TRISpH(double T, double P)
{
	GeoSystem system = AssembleSystemTRIS();
	
	GeoState state(system);
	
	state.SetTemperature(T);
	state.SetPressure(P);
	
	GeoEquilibrator equilibrator(system);
	
	cout << setw(15) << left << "mTRIS";
	cout << setw(15) << left << "mNaCl";
	cout << setw(15) << left << "pH";
	cout << setw(15) << left << "ln(gTRIS)";
	cout << endl;

	BOOST_FOREACH(double mTRIS, Range(1.0E-6, 0.02, 0.005)) 
	{
		BOOST_FOREACH(double mNaCl, Range(1.0E-6, 5.0, 1.0))
		{
			equilibrator.ResetConditions();
			
			equilibrator.ImposeChargeBalance();
			equilibrator.ImposeSpeciesMass("H2O", 55.508);
			equilibrator.ImposeElementMass("Na", mNaCl);
			equilibrator.ImposeElementMass("Cl", mNaCl);
			equilibrator.ImposeElementMass("TRIS", mTRIS);
		
			equilibrator.Equilibrate(state);
			
			cout << setw(15) << left << mTRIS;
			cout << setw(15) << left << mNaCl;
			cout << setw(15) << left << state.Get("pH");
			cout << setw(15) << left << std::log(state.Get("ActivityCoefficient", "TRIS"));
			cout << endl;
		}
		
		cout << endl;
	}
}

void TRISpHTRange(const vector<double>& rangeT, double P, double mNaCl, double mTRIS)
{
	GeoSystem system = AssembleSystemTRIS();
	
	GeoState state(system);
	
	state.SetPressure(P);
	
	GeoEquilibrator equilibrator(system);
	
	cout << setw(15) << left << "T";
	cout << setw(15) << left << "pH";
	cout << setw(15) << left << "ln(gTRIS)";
	cout << endl;

	cout << fixed << setprecision(5);
	
	BOOST_FOREACH(double T, rangeT) 
	{
		state.SetTemperature(T);
		
		equilibrator.ResetConditions();

		equilibrator.ImposeChargeBalance();
		equilibrator.ImposeSpeciesMass("H2O", 55.508);
		equilibrator.ImposeElementMass("Na", mNaCl);
		equilibrator.ImposeElementMass("Cl", mNaCl);
		equilibrator.ImposeElementMass("TRIS", mTRIS);

		equilibrator.Equilibrate(state);

		cout << setw(15) << left << T;
		cout << setw(15) << left << state.Get("pH");
		cout << setw(15) << left << std::log(state.Get("ActivityCoefficient", "TRIS"));
		cout << endl;
	}
}

//const Multiphase AssembleTRISMultiphase()
//{
//	vector<string> aqueousSpecies = {"H2O", "Na[+]", "Cl[-]", "H[+]", "TRIS", "OH[-]", "NaCl(aq)", "HCl(aq)", "NaOH(aq)", "TRISH[+]"};
//	
//	ElectrolyteSolution sol(aqueousSpecies);
//	
//	AqueousActivity aqueousActivity(aqueousSpecies);
//	
//	aqueousActivity.SetActivityCoefficient("H[+]", PitzerHydronModel(sol));
//	aqueousActivity.SetActivityCoefficient("TRIS", PitzerTRISModel(sol));
//	aqueousActivity.SetActivityCoefficient("TRISH[+]", PitzerTRISHModel(sol));
//	
//	Phase aqueousPhase;
//	
//	aqueousPhase.SetName("Aqueous");
//	aqueousPhase.SetSpecies(aqueousSpecies);
//	aqueousPhase.SetActivity(aqueousActivity);
//	
//	Multiphase multiphase({aqueousPhase});
//	
//	return multiphase;
//}
//
//const ReactionSystem AssembleTRISReactionSystem(const Multiphase& multiphase)
//{
//	EquilibriumConstant K(EqulibriumConstantTRISH);
//	
//	ReactionEquation equationTRISH({"TRISH[+]", "TRIS", "H[+]"}, {-1.0, 1.0, 1.0});
//	
//	Reaction reactionTRISH(equationTRISH, K, multiphase);
//	
//	vector<Reaction> reactions;
//	
//	reactions.push_back(AssembleReaction("OH[-]", multiphase));
//	reactions.push_back(AssembleReaction("NaCl(aq)", multiphase));
//	reactions.push_back(AssembleReaction("HCl(aq)", multiphase));
//	reactions.push_back(AssembleReaction("NaOH(aq)", multiphase));
//	reactions.push_back(reactionTRISH);
//	
//	ReactionSystem eReactions(reactions);
//	
//	return eReactions;
//}

const Multiphase AssembleTRISMultiphase()
{
	vector<string> aqueousSpecies = {"H2O", "Na[+]", "Cl[-]", "H[+]", "TRIS", "OH[-]", "TRISH[+]"};
	
	ElectrolyteSolution sol(aqueousSpecies);
	
	AqueousActivity aqueousActivity(aqueousSpecies);
	
	aqueousActivity.SetActivityCoefficient("H[+]", PitzerHydronModel(sol));
	aqueousActivity.SetActivityCoefficient("TRIS", PitzerTRISModel(sol));
	aqueousActivity.SetActivityCoefficient("TRISH[+]", PitzerTRISHModel(sol));
	
	Phase aqueousPhase;
	
	aqueousPhase.SetName("Aqueous");
	aqueousPhase.SetSpecies(aqueousSpecies);
	aqueousPhase.SetActivity(aqueousActivity);
	
	Multiphase multiphase({aqueousPhase});
	
	return multiphase;
}

const ReactionSystem AssembleTRISReactionSystem(const Multiphase& multiphase)
{
	EquilibriumConstant K(EqulibriumConstantTRISH);
	
	ReactionEquation equationTRISH({"TRISH[+]", "TRIS", "H[+]"}, {-1.0, 1.0, 1.0});
	
	Reaction reactionTRISH(equationTRISH, K, multiphase);
	
	vector<Reaction> reactions;
	
	reactions.push_back(AssembleReaction("OH[-]", multiphase));
	reactions.push_back(reactionTRISH);
	
	ReactionSystem eReactions(reactions);
	
	return eReactions;
}

const GeoSystem AssembleSystemTRIS()
{
	Multiphase multiphase = AssembleTRISMultiphase();
	
	ReactionSystem eReactions = AssembleTRISReactionSystem(multiphase);
	
	GeoSystem system(multiphase, eReactions);
	
	return system;
}
